Complex Network Analysis of State Spaces for Random Boolean Networks 
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We apply complex network analysis to the state spaces of random Boolean networks (RBNs). 
An RBN contains A*" Boolean elements each with K inputs. A directed state space network (SSN) 
is constructed by linking each dynamical state, represented as a node, to its temporal successor. 
We study the heterogeneity of an SSN at both local and global scales, as well as sample-to-sample 
fluctuations within an ensemble of SSNs. We use in-degrees of nodes as a local topological measure, 
and the path diversity ^ of an SSN as a global topological measure. RBNs with 2 < K < 5 
exhibit non-trivial fluctuations at both local and global scales, while K = 2 exhibits the largest 
sample-to-sample, possibly non-self-averaging, fluctuations. We interpret the observed "multi scale" 
fluctuations in the SSNs as indicative of the criticality and complexity of K = 2 RBNs. "Garden 
of Eden" (GoE) states are nodes on an SSN that have in-degree zero. While in-degrees of non-GoE 
nodes for K > 1 SSNs can assume any integer value between and 2^, for K = 1 all the non-GoE 
nodes in an SSN have the same in-degree which is always a power of two. 

PACS numbers: 05.45.-a, 89.75.-k, 89.75.Fb, 89.75.Da 



I. INTRODUCTION 

In this paper we apply complex network analysis to dis- 
crete, disordered, deterministic dynamical systems. The 
set of all trajectories of such a system can be described as 
a directed network. Each dynamical state is represented 
by a node and linked to its unique temporal successor 
by a directed link, giving the state space network (SSN). 
Thus the out-degree of a node is one. The irreversibility 
of the dynamics implies that a node can potentially have 
many in-coming links. The number of in-coming links at 
a node, or its in-degree, can vary from node to node, de- 
pending on the dynamical system considered. A wide 
dispersity of degrees characterizes many complex net- 
works ^2]. Therefore, complex network analysis may offer 
a useful alternative to traditional analyses of dynamical 
systems, such as the characterization of spatiotemporal 
patterns [l,l3,S[EI3l- The first results of such an analysis 
on one dimensional cellular automata (CA) showed that 
heterogeneity of the SSNs at both the local and global 
scales distinguishes "complex" dynamics from simple dy- 
namics 1]. Here we exploit complex network theory to 
characterize disordered dynamical systems by examining 
SSNs for ensembles of random Boolean networks. 

Random Boolean networks (RBNs) were introduced by 
KaufFman ^ as models of gene regulation and have been 
extensively studied over the years by physicists as exam- 
ples of strongly disordered systems [1, [13, [HI, [H, [ll]- 
The dynamics of each of the N Boolean elements in an 
RBN is given by a Boolean function of K randomly cho- 
sen input elements. Different realizations of the input 
elements and the Boolean functions for an RBN lead to 
an ensemble of RBNs for a given {N, K). 

In the thermodynamic {N oo) limit RBNs exhibit a 
phase transition between chaotic and frozen phases pass- 
ing through a critical phase for K = 2 [3| . In the frozen 



phase, K < 2, the Hamming distance between two per- 
turbations of the same state quickly die out. On the 
other hand, in the chaotic phase, K > 2, perturbations 
grow exponentially in time. Substantial analytical work 
has focused on the number of attractors as well as their 
lengths fisl . [l6l . [iTj . which have been found exactly for 
K = 1 |l8l |. Krawitz and Shmulevich computed the 
entropy of basin sizes and showed that for critical RBNs 
it scales with the system size. Otherwise, it asymptotes 
to a constant. 

In unrelated developments, a variety of statistical 
methods have been established to characterize the struc- 
ture of complex networks. Measures include probabil- 
ity distributions for the node degree, clustering (the ten- 
dency of nodes that share a common neighbor to also 
be linked to each other directly), motifs (overrepresented 
subgraphs in the network), etc. d, [13, [Hj- Many real 
world networks such as regulatory networks [22|, the 
world-wide web [23| , or the correlation structure of earth- 
quakes [l^, [2^ ) differ markedly from a random graph - 
where the degree distribution is Poisson and clustering is 
absent. They often display "fat-tailed" or even scale- free 
degree distributions. 

In this paper we show that strong deviations from ran- 
dom graph behavior also occur in the SSNs of RBNs. In 
particular, for sufficiently small K, the probability dis- 
tribution for the number of incoming links to nodes in 
the SSN, P{k), has a fat-tail. In contrast, a randomly 
rewired null model for the SSN (i.e., a random map) has 
Poissonian in-degree distributions. Also for small K, the 
maximum in-degree of any node in the SSN, kmax, dis- 
plays scaling behavior with respect to size of the SSN, 
M ~ 2^ . While kraax is a measure of local heterogene- 
ity, we use thepath diversity, T), to characterize global 
heterogeneity [l| . We show that T> grows linearly with N 
for K = 1, while it grows faster than linearly with N for 
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2 < -ff < 5. For if > 6, r> scales with SSN size M. In 
addition, for K = 2 (where RBNs are critical) the sample 
to sample fluctuations of V are the largest and might be 
non-self-averaging [2^, [23|. We speculate that SSN fluc- 
tuations at these three different scales (local, global, and 
sample to sample) for X = 2 RBNs are associated with 
criticality in the thermodynamic limit M oo. 



A. Summary 

In Sectioniniwe discuss the procedure used to construct 
RBNs, and their SSNs. The in-degree and the path diver- 
sity are also defined. In Section UlIAl we present results 
for the behavior of the nodes' in-degrees in an ensemble 
of SSN. In SectionHnXll we discuss the SSNs for if = 1 
RBNs, which have unique features not shared by other 
if s. In Section IIII Bl we examine the behavior of the 
path diversity. Discussion of our results and concluding 
remarks are found in Section IIVI 



II. DEFINITIONS 

An RBN consists of N Boolean variables ct^ G 0, 1 with 
i = 1 ■ ■ ■ , N . The dynamics of each element is determined 
by a Boolean function of K randomly chosen input ele- 
ments, 

(T,{t + l) = h{<j,,{t),<j,^{t),...,a,^{t)) , (1) 

where ai - {t) is the value of the input to (Ji at time t. 
The function fi is randomly chosen to be 1 with proba- 
bility p and with probability 1 — p for each set of values 
of its arguments. We consider only the case of unbiased 
RBNs with p — 1/2, except when K = 1, where we study 
two different cases. For if = 1 there are four possible 
boolean functions for each element. Instead of choosing 
them uniformly, using only the copy (/(ct) = a) or the 
invert function (/(cr) — NOT(cr)) for each a i (with equal 
probability), leads to a critical if = 1 RBN Q. 

We use synchronous update for the dynamics, i.e. all 
the Boolean elements in the network are updated in par- 
allel at each time step. The networks are set up by 
choosing if different random inputs for each ai. While 
allowing self-connections, we do not allow multiple con- 
nections. We do not impose any connectivity constraint 
on the RBNs. 

An RBN with N elements has A/" = 2^ different dy- 
namical states. These states are nodes of a directed net- 
work. A link from A to B indicates that A evolves to 
B in one time step, making B the image of A, or A a 
pre-image of B. This directed network forms the state 
space network (SSN). The SSN typically consists of dis- 
connected clusters, or basins of attraction. Each basin 
contains transient states, which are visited no more than 
once on any dynamical trajectory, and attractor states 
that may be visited infinitely often. Garden of Eden 



states (GoE) are transient states that cannot be reached 
from any other state, i.e they have no pre-images. Ex- 
amples of some state space clusters are shown in Fig. [T] 
for = 9 and different values of if. 




FIG. 1: (Color online). One connected cluster of an SSN for 
RBNs with A'' = 9 and different values of if displayed using 
the program 'Pajek' [28l ]. Nodes on the attractors are drawn 
in red; the other colors indicate distance from the attractor. 
For instance, for if = 1, green nodes are distance one from 
the attractor, blue nodes are distance two and magenta nodes 
are distance three. Note that for if = 1 all nodes are either 
Garden of Eden states or hubs with all the hubs having the 
the same in-degree. RM stands for a random map. 

Random maps are the limit of RBNs when N,K ^ 
oo 29]. By construction, the state space of a random 
map forms an Erdos-Reyni random graph with a Poisson 
in-degree distribution, P{k) = e~^/k\ with mean (k) = 
1. We constructed random map SSNs by picking the 
image of each node uniformly randomly from the JV = 2^ 
possible nodes. 

Identifying all attractors states and lumping them into 
a single node turns the SSN into a rooted tree. To study 
the global heterogeneity of this tree, we use path diversity 
v. It quantifies topological variations in the paths con- 
necting the GoE states to the root. It is similar to other 
global measures of tree like structures such as tree diver- 
sity [3^ and topological depth Q . V measures the num- 
ber of non-equivalent choices encountered when following 
each reversed path from the root to the GoE states. 

Specifically, V is computed as follows: each GoE state 
is assigned a diversity equal to one; a state with a single 
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pre-image has the same diversity as its unique pre-image; 
and the diversity of a node with more than one incoming 
hnk is the sum of all distinct diversities of its pre-images 
plus one. For instance, if a node has in-degree 6, and 
three of its pre-images have diversity 4, two have diversity 
5, and the last has diversity 8, then that node's diversity 
is 4 + 5 + 8 + 1 = 18. Finally, the path diversity V of the 
entire SSN is the diversity of the root. 

Table U provides a summary of the notation used. Note 
that we refer to nodes of the RBNs as "elements" and 
otherwise reserve the phrase "node" only for SSNs. 



RBN 


random Boolean network 


iV 


number of elements of an RBN 


K 


number of inputs to each element of an RBN 


SSN 


state space network of an RBN 


JV 


number of nodes in an SSN (TV = 2^^) 


k 


in-degree of a node in an SSN 


kmax 


largest in-degree of a node in an SSN 


V 


path diversity 



TABLE I: Notations used in this paper. 



III. RESULTS 
A. In-Degree 

An elementary description of local heterogeneity of 
a network is its degree distribution. We computed in- 
degree distributions, P{k), of the SSNs. These SSNs were 
obtained for RBNs with 1 < < 10 and 10 < iV < 24, 
as well as random maps of system sizes 10 < < 24. For 
iiT = 1 we distinguish between RBNs constructed using 
all four boolean functions and the critical K = 1 RBN 
discussed earlier. In Section IIII A II we discuss the in- 
degrees ioi K — 1 RBNs, where we obtain exact analytic 
results. We present in Section IIII A 21 numerical results 
for > 2 RBNs. 



1. K—1 Networks 

We start with the case oi K = 1 critical RBNs, where 
all functions fi{(Ji-^) are either copy or invert. If none of 
the boolean elements is a leaf on the RBN, then the value 
of every element ai at the next time step is determined 
by the current value of of some element aj . If the state 
S differs from S' in the value of cr^, their pre-images will 
differ in the value of aj. Thus on a = 1 critical RBN 
without leaves, the mappings are one-to-one, there are 
no GoE states and each hub has in-degree one. 

When the RBN contains L leaves in addition to Q non- 
leaf elements, the state of the L + Q elements can be 



written as a concatenation, 

S = L©Q. (2) 

If ai is among the L leaves, it cannot effect the next 
step value of any element. The next value of itself is 
determined by the current value of a non-leaf element aj . 
The image of S under the RBN rules will then be of the 

form, 

Z(S) = .F(Q) ® g(Q), (3) 

where gives the next-step state of the leaves, and is 
a function only of the non-leaves. This means that at 
least 2^ states that differ only in the values of the L leaf 
elements will be mapped to the same image, i.e. the in- 
degree of each non-GoE state ("hubs"), is kh > 2^. To 
see that kh = 2^ , we note that 

X(S') ^ Z(S) if Q' ^ Q. (4) 

Consider two states Q and Q' of the non-leaf elements 
that differ in the state of the element a a- If ctq is an 
input to at least one non-leaf element, then t/(Q') ^ 
fJ(Q). If (Ta is an input to at least one leaf element, 
then 3-{Q,') ^ •?^(Q)- Since CTq is not a leaf, at least 
one of these two cases has to be true, proving Eq. Q. 
Thus the cardinality of the set {Z(S)} is at least equal 
to the cardinality of the set {Q}, which is 2^ . If for 
each hub kh > 2^, the sum of hub in-degrees > 2*^2^ = 
2Q+L ^ 2^. Since there are only 2^ nodes in the SSN, 
this inequality has to be saturated, giving us that kh = 
2^. 

The boolean function used for K = 1 critical RBNs 
are either copy or invert. For the general K — 1 case, 
there are four possible functions, two of which are copy 
and invert. The other two are constant functions, which 
map all aruguements to one value (either zero or one). 
If an element at is input only to elements with constant 
functions, it behaves effectively like a leaf. We will refer 
to such elements along with the leafs, as effective leaves. 
Using this observation we can extend the above argument 
for if = 1 critical RBNs to the general case of if = 1 
RBNs. Thus all hubs in the SSN of a if = 1 RBN have 
the same in-degree of the form kh — 2', where I is the 
number of effective leaves. However, kh will vary over 
the different realizations of the RBN. 

Construction of an instance of a if = 1 RBN can be 
considered as N rolls of an iV-faced die. Each face rep- 
resents an element. The outcome j of the r-th roll is the 
input to the r-th element of the RBN. The function f-r 
in Eq. ([T]) that determines the updates of the r-th ele- 
ment is chosen randomly over the function's arguments, 
aj in the present case. Since for if = 1 there is only 
one input, fr will be a constant function if it maps both 
values of the input (0,1) to the same value, either 1 (with 
probability p^) or (with probability (1 — p)^). There- 
fore, the probability that the chosen function is not a 
constant function is q — 1 — — (1 — p)^. We flip a coin. 
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after each roll, to determine if the function fr is a con- 
stant function or not. If it is not a constant function, we 
mark the displayed face of the die, and keep track of the 
number of marked faces until t = N. The marked faces 
exclude the candidates for effective leaves. After N die 
rolls and coin flips, the number of marked faces will be 
m = N — I, where I is the the number of effective leaves 
in the constructed RBN. The distribution of m evolves 
from the r-th to the r + 1-th roll according to. 



Mr + l,m = 9 1 



1 



N 



Mr,m-l+(l-q{l 



^)) M..™, 
(5) 



under the boundary conditions that 

Mr^m = 0, if TO < or m > t, 

and the initial condition Afo,o = 1- 

Thus q = 1 in Eq. ([5]) generates the distribution of 
leaves for K ~ 1 critical RBNs. In this case, the solution 
to Eq. ^ is 



Mw.m - 77F(^)'^^^"'"^■' 



(6) 



where CN,m are Stirling numbers of the second kind, and 
enumerate partitions of an iV-set into m non-empty sub- 
sets ^H. For q ^ 1 analytical expressions are hard to 
obtain. However, solving Eq. ([5]) numerically is straight- 
forward. 

We can use the distribution of to to generate the distri- 
bution of logjfc/i — N — m. In Figs. [2] and [3] we present 
the distribution of logarithm of hub in-degrees logjfc/i for 
the K = 1 critical and K = 1 non-critical RBNs respec- 
tively. The hub-degrees were obtained for an ensemble 
of 2, 000 randomly chosen RBNs with 10 < TV < 24. We 
also show our analytical results obtained from Eq. ([5]) 
with g = 1 for the K = 1 critical RBNs and with q — 1/2 
for the fs: = 1 RBNs. 

We can use Eq. ([5]) to obtain the mean number of ef- 
fective leaves, 



(to). 



mMr+l,m 

T+1 

E 



m=0 

T+1 

m=0 

which reduces to, 

{m)r+i = q - 



TO — 1 
9 11 — I Mr.,n-1 



1 - (7(1 - -) ) M,,, 



N 



{m)r- 



This recursion gives the solution for the mean number 
of effective leaves N — {m)N, which is also the mean of 
log2^/i, 



(log^fc/.) = TV (l - I) 



N 



Ne 



as N oo. (7) 




FIG. 2: (Color online). The PDF P(log2 k) of the log of the 
hub in-degree, logj kh for K = 1 critical RBNs. The three 
dashed lines with empty symbols were obtained using Eq. (|5| 
with g = 1. 
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FIG. 3: (Color online). The PDF P(log2 k) of the log of the 
hub in-degree, logj kh for K = 1 RBNs. The three dashed 
lines with empty symbols were obtained using Eq. ((5} with 
q = 1/2. 



2. Networks with K >2 

For K > 1 RBNs the boolean functions are not the 
simple copy /invert and constant functions as in the case 
oi K — 1. For example, consider the truth table in Table 
im (Ti and a2 are inputs to 0-3 in a if = 2 RBN. (T2 effects 
the value of only when ai — I. Conversely, <ti effects 
the value of CT3 only when (72 = 1- Such interactions 
between the inputs ior K > 1 RBNs make it harder to 
extend the analytical aruguements that we applied to the 
K — 1 case. 

Fig. H presents P(fc) for SSNs with iV = 18 (TV = 
2^^). It is broad for K = 1, becomes narrower with 
increasing K, and converges for if ^ A'' to the in-degree 
distribution of the random map. The random map in- 
degree distribution itself converges to a Poissonian with 
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TABLE II: A function iox K = 2 RBNs 



{k) = 1 when N ^ oo. Note that k for if = 1 is either 
zero for the GoE states, or takes the same power of two 
value for the hubs. 
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FIG. 4: (Color online). The in-degree distribution function, 
P{k), for various connectivities, K, of the RBN, on a log-log 
scale. The distribution is computed for RBNs of size A'^ = 18 
by averaging the PDFs over 200 different realizations for each 
K. We used fc -I- 1 instead of k on the a;-axis, so that we 
could also show values for fc = 0. The black solid line shows 
a Poissonian with (fc) — 1 for comparison. K = IC refers to 
K = 1 critical RBNs. 



The in-degree distribution P{k) does not exhibit eas- 
ily described behavior such as scahng. However, the re- 
sults of Ref. [H suggest that a more useful quantity is 
the largest in-degree kmax- Scaling behavior in kmax was 
shown to be a necessary but not sufficient signature of 
complex dynamics. Fig. [5^ shows that 



(log. 



2 ^max 

) i^rN for K <6, 



(8) 



where the angular brackets indicate an average over dif- 
ferent realizations of the RBN. The exponents vk appear 
to obey the relation: 



i^K = -0.07K + 0.68 for 1 < K < 6, 



(9) 



where the error is ±1 in the last digit. For the system 
sizes studied, it is not possible to determine the asymp- 
totic limit for larger K. However, the behavior for large 
K approaches the random map result. 

For the random map, P(fc) tends to a Poisson distribu- 
tion when N ^ oo. If A/" different values of k are chosen 



from a distribution P{k), the expected 
can be related to J\f by 



1 

77- 



(10) 



For a Poissonian with mean one, 

kl ~ N" 



E 



This gives the following bounds on k„ 



< 



< 



E 



m—0 ^ 



The sum in the last term in the equation evaluates to 
1 — 1/kmax which goes to one for large kmax, allowing us 
to write 



1 



Taking logarithms and using Stirling's approximation, 
the above equations yields for large M 



. In k„ 



lnA/' = iVln2. 



Finally we take logarithms again to get, 

log2 kmax ~ loga N. 



(11) 



(12) 



As shown in from Fig. [5^, (log2 kmax) is very well de- 
scribed by Eq. (fT^ for random maps and RBNs with large 
K in the limit N oo. 

In order to see whether kmax/N is self averaging, 
i.e. whether the fluctuations of kmax I N become neg- 
ligible for large N , we also plot in Fig. O the ratio 
A^^^ log2(fcmaa:)- A first look at Fig. [5] might suggest 
that both ways of averaging lead indeed to the same re- 
sults, and self-averaging is satisfied. That this is not the 
case is demonstrated in Fig. [SI where we show the ratio 
(log2 kmax) I log2 {kmax) vcrsus A^. Wc scc from this figure 
that {kmax) scales with AC as {kmax) ~ A/"^^ for if < 6, 
i.e. log2{kmax) ^ trN, similarly to Eq.ijS]). However, 
tk 7^ i^K- For large K it seems that that the fluctua- 
tions are less important, (logj fcmaa;)/ log2 (fcmaa;) 1 for 
K oo. 

The sample-to-sample fluctuations in kmax a-i'e cap- 
tured by its probability distribution, shown in Fig. [T) 
We plot the probability distribution of log kmax multi- 
plied by its standard deviation, a(\og2 kmax)- The mean 
and standard deviation were computed independently for 
each A^ and K. The plots are well-fitted by a Gaussian 
for 2 < AT < 6, which suggests that P{kmax) is a log- 
normal, for sufficiently large A^. For AT > 6, the plots 
indicate deviations from log-normal behavior. 
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FIG. 5: (Color online), (a) The ratio 



(log 



is plotted 

against A'^. (b) The ratio '"^^^^"""'^ is plotted against A''. 
Both figures show systematic deviations from scaling for K > 
6. The solid lines in (a) and (b) are '"^"'^ and were plotted 
for comparison. The statistics, in this figure and Fig. |6] were 
obtained by sampling 2000 different realizations for each A'' 
and K. 



B. Path Diversity 

Scaling of kmax with system size for K < 6, indicates 
asymptotic local heterogeneity of the SSNs. As shown in 
Ref. Il| local heterogeneity is often insufficient to distin- 
guish simple from complex dynamics. A complimentary 
global topological measure, the path diversity, detects 
global heterogeneity in SSNs. In Ref. [l| it was found 
that complex cellular automata show scaling behavior in 
both local and global measures. Here we investigate the 
path diversity of SSNs of RBNs. 

Fig. [5] suggests that (logj V) ^ CkN for K > 5 and for 
the random map. However, the plots are not conclusive 
as to whether or not (log2 T>) becomes linear in N for 



FIG. 6: (Color online). The ratio between the mean of 
the log and the log of the mean of the largest in-degree, 
y = io7/(fc"°^) . is plotted against the size of the RBN, N. 
If the two quantities were scaling with the same exponent, 
all curves should tend to y — 1 for large A''. The plot shows 
that (log2 kmax) and \og2{kmax) are scaling with different ex- 
ponents for K <Q. 



K < 5, for large system sizes, or with large finite size 
corrections. Fig. [5] shows Aiog^ jv (logj 2?) , the discrete 
derivative of (logj T)) with respect to log2 A''. For K = 1, 
the figure indicates that (log2 T>) logj A^. A faster than 
logarithmic growth of (log2 T>) is seen for AT = 2. 

An example of a simple tree that exhibits logarithmic 
behavior of path diversity is the Cayley tree. On such 
an SSN each transient state has exactly z pre-images, 
except for the root which has z + 1 pre-images. Given the 
symmetry of the Cayley tree, all the nodes at the same 
distance from the root will have the same path diversity. 
Furthermore, the path diversity is incremented by one 
with each step towards the root. This makes the path 
diversity equal to the depth of the tree which can be 
expressed in terms of the number of nodes JV: 



V = \og,{{z-l)M+l}. 



(13) 



For a maximally heterogeneous SSN each transient 
node will contribute to the path diversity. The path di- 
versity of an SSN is the diversity of the root. Diversity 
of a node is the sum of all distinct diversities of its pre- 
images plus one if it has more than one pre-image. Thus 
the diversity of a node is bounded above by, 

j-<i 

where j < i indicates that j is a pre-image of i. In the 
general case, only distinct values of the diversities will 
contribute to the sum on the right hand side. We can 
iteratively use this upper bound to calculate the path 
diversity of an SSN. Denoting the root of the SSN by •, 



V < I 
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(a) 




FIG. 7: (Color online). The rescaled PDFs of the log of the 
largest in-degree P(log2 kmax) for {a,)K — 2, 4, and 6, and 
(b) K — 7, 9, and the random map. The dashed lines are 
Gaussian distributions with mean zero and variance one. The 
plots indicate that kmax is distributed according to a log nor- 
mal distribution for K = 2, 4 and 6. Similar results hold for 
K — 3 and 5. Deviations from the Gaussian distribution are 
seen for K > 6. In (a) the distributions for K = 4 and K = 6 
were offset by 7.5 and 15 units on the x-axis for clarity of pre- 
sentation. In (b) the distributions for K — 9 and the random 
map were offset by 10 and 20 units. 



= 1 + Ei+E E^»= 

ii-<» 12^*1 13^12 i2^il il'i* 

and so on. The subscript h in ih represents the distance 
of the node ih from the root. If ih is a GoE state, its 
diversity is defined to be one. The first term on the right 
hand side of the last line counts the root, the first sum 
counts the root's pre-images, and the two nested sums 
count the pre-images of the root's pre-images. Iteratively 



FIG. 8: (Color online). Log path diversity as a function of 
the system size. {log2'D) /N , is plotted as a function of A'', 
for various values of K and the random map. For K > 5, 
(logj 2?) asymptotes to CkN, while the plots have the opposite 
curvatures for K < 4. 
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log^N 

FIG. 9: (Color online). The discrete derivative of (logj ©) 
with respect to log2 A'', Aiog^ jv (log2 f ) , is plotted against 
logj A'^ for K — 1,2 and K = 1 critical. Aiog2 jv (logj O) is 
defined as ^'"sa Pa) — (iog2 p.-i) ^ ^ horizontal line of this plot 

indicates a relationship of the form (logj V) ~ logj A'^. For the 
system sizes studied, the data suggests a logarithmic growth 
of (logj 23) with respect to A''. For K — 2, Aiog^ jv(log2 23) 
show growth trends which are indicative of a faster than log- 
arithmic growth of (log2 23) with A'. The data is inconclusive 
for K = 1 critical. 



continuing this calculation we can see that the right hand 
side evaluates to M , the size of the SSN. Thus we have 
as the upper bound to the path-diversity, 

V<N^2^, (14) 

or that, 

\0g2V < N. (15) 
This bound is nearly reached for binary tress where 
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single leaves branch off at every branching point. 

The logarithmic growth of {log2 T>) for if = 1 net- 
works suggests a different class of dynamical complexity 
than for K >2 networks. Meanwhile, random map-like- 
scaling of (log2 kmax) for RBNs with sufficiently large K 
suggests relatively simple dynamics. However, on the ba- 
sis of those two measures, one cannot cleanly distinguish 
the behavior oi K — 2 critical RBN from that oi K > 2 
RBNs. We now show that sample-to-sample fluctuations 
in log2 V enable such a distinction. 
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FIG. 11: (Color online). The variance of the log of the path 
diversity ^^(logjl') is plotted as function of the size of the 
RBN, N, for various values of K. cr^ (loggia) shows non- 
monotonic behavior as a function of K, it grows the fastest 
for K — 2. On the other hand, (T^(log2 D) appears to tend to 
a constant or decreases with increasing A'^, for large values of 
K. For the random map, (T^(log2 23) is a decreasing function 

of A'^. The inset shows a plot of Ti{T>) = "(^^"^^1 as function 
of A'^. The data for K — 2 suggests that logj 23 might be 
non-self- averaging in the thermodynamic limit. 



FIG. 10: (Color online). The log of the PDF of the log of 
the path diversity, P(log2 23), for various values of K and 
N = 24. P(log2 V) is the broadest for K ^ 2. P(log2 V) 
becomes narrow for large values of K. 



We report the distribution of log2 1?, P{\og2'D), in 
Fig. [TUl For a fixed system size, P(log2 T>) is broadest 
for if = 2 RBNs. To quantify the width of P(log2l>), 
we study its variance. Fig. [TT] shows that the variance 
of log2 23, cr^ (log2 23), grows fastest with N for K = 2. 
Unlike {log2 kmax) and (loggP), O'^(log2 23) shows non- 
monotonic behavior as a function of K . We interpret this 
as an indication that the criticality of if = 2 RBNs can 
be associated with large sample-to-sample fluctuations of 
its SSNs. 

Generally, large sample-to-sample fluctuations in dis- 
ordered media may lead to the absence of self-averaging. 
A canonical measure of self-averaging is the ratio: 

A system is said to be self-averagin g w ith respect to 23 
if 7?.(23) goes to zero for large N |26l. |27|: otherwise, it is 
said to lack self-averaging. For the system sizes we are 
able to study, the evidence for or against self-averaging 
is not completely conclusive. Nevertheless, as indicated 
in the inset of Fig [Til K = 2 RBNs show the largest 
values for TZ and are therefore the most likely to exhibit 
non-self-averaging behavior in the thermodynamic limit 
of large system size. 



IV. DISCUSSION AND CONCLUSION 

We have studied the topology of state space networks 
(SSNs) for ensembles of random Boolean networks. Each 
dynamical state of the Boolean network corresponds to 
a node in the SSN and is linked to its successor state. 
We characterize the heterogeneity of these SSNs at the 
local, node, scale by the distribution of the in-degrees 
and the scaling of the largest in-degree. Global hetero- 
geneity over all paths in the SSN is characterized by the 
path diversity. For elementary 1-d cellular automata, 
it was demonstrated in Ref. tJ that simultaneous scal- 
ing behavior in both kmax and 23 indicates "complex" 
spatio-temporal dynamics (Wolfram class IV and partly 
class III). On the other hand, it was found in fT] that one 
or both of these does not scale for CA in class I or II - 
which all have simple dynamics. 

As is well-known, RBNs exhibit a phase transition 
between chaotic (if > 2) and frozen (if < 2) behav- 
ior, if = 2 RBNs are critical and therefore more com- 
plex than RBNs with K ^ 2 which mainly show frozen 
{K < 2) or chaotic (if > 2) behavior [32]. In fact, RBNs 
in the frozen phase (e.g. if = 1) resemble class II CA as 
they almost always rapidly go to one of many attractors 
with a short period. In addition, RBNs in the chaotic 
phase resemble some class III CA in that they have long 
transients and attractors with large periods. We have 
investigated whether or not the different phases of RBNs 
can be distinguished on the basis of a topological analysis 
of the corresponding ensembles of SSNs, and the extent 
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to which these phases overlap in the behavior of their 
SSNs. 

Man y a nalytical results are known for K = 1 
RBNs [il mi (and also for the random map, ^29^). In 
the present paper we derived analytical results for the 
structure of the state space network for K — 1 RBNs, 
and verified them using numerical methods. We related 
the hub in-degree kh to I, the number of effective leaves 
in the RBN. Simple arguments show that all the hubs 
have the same in-degree fc/j = 2K We also presented a 
stochastic process in Eq. (O that generates the distribu- 
tion of I for any if = 1 RBN. For = 1 critical RBNs 
we could obtain an analytical expression for the distri- 
bution of the hub in-degree kh , and compared numerical 
solutions to Eq. ^ with the simulations for the general 
K = 1 RBNs. 

Using Eq. (O we derived that ior K ^ 1, logjfc/i scales 
with the size of the RBNs as (logjfc/j) = e~'^N, where q 
is the fraction of non-constant functions. Numerically we 
find that for 2 < A' < 6, log2 kmax also exhibits scaling 
as a function of the size of the RBN, (logj k^ax) ~ vj^N . 
The scaling exponent vk is largest ioi K = 1 and mono- 
tonically decreases for larger values of K. For K = 1 the 
hub in-degree kh depends exponentially on the number 
of effective leaves, I, which scales linearly with N . The 
state of effective leaves is forgotten after the first time 
step during evolution of the RBN. For if > 2, there are 
more than one inputs to each element. Whether the state 
of an element is forgotten or not depends on the state of 
other elements in the RBN. Observed linear scaling of 
\-Og2kmax with N for 2 < A'^ < 6 indicates an explanation 
similar to the one we presented for K = \. However, more 
analytical work is required to fully understand the rela- 
tionship between k^ax and the structure of the RBN for 
K > \. For larger AT, we notice that the K,N ^ 00 limit 
of the RBNs is the random map [1^ . We have shown an- 
alytically that for the random map, log2 kmax ^ log2 N . 

On the other hand, our numerical results indicate that 
(log2 V) scales with the size of the SSN only for suffi- 
ciently large values oi K {K >^), (log2 V) ~ C,kN. C,k is 
largest for the random map and monotonically decreases 
with decreasing K (which is the opposite oIvk)- The log- 
arithmic scaling with N of (log2 T)) ioi K — 1 translates 
to a logarithmic scaling of V with the size of the SSN, 



TV. This compares with the path diversity of a Cayley 
tree in Eq. p^ . On an SSN with Cayley tree structure, 
all nodes, except GoE nodes, have the same in-degree. 
In fact we have shown the same for if = 1 SSNs. For 
2 < if < 5, it is not clear from the data whether (log2 V) 
scales linearly with N . It is clear, however, that (log2 T>) 
grows faster-than-linear with log2 N . 

Together, the absence of scaling for k,nax and V cor- 
rectly rule out the random map (and most likely RBNs 
with large K) and K = 1 RBNs, from reporting high 
dynamical complexity. However, neither of these mea- 
sures addresses sample-to-sample fluctuations, which are 
generally important in the characterization of disordered 
systems. We find that the probability distribution func- 
tion of kmax converges to a log-normal distribution for 
2 < if < 6. While this is also seen for log2 V, unlike the 
monotonic dependence of the variance of log2 k^ax , the 
variance of log2 V does not vary monotonically with if. 
In fact, the variance grows the fastest with the size of the 
SSN for critical if = 2 RBNs. Numerical results also sug- 
gest that log2 V may possibly exhibit non-self-averaging 
behavior for if = 2. 

Our results support the conclusion that heterogeneity 
in SSNs can be associated with complex dynamics in dis- 
ordered systems, if = 2 RBNs are distinguished from 
other RBNs by simultaneously exhibiting three kinds of 
network heterogeneity. The first is a local heterogeneity 
on the node level indicated by the scaling of log2 kmax 
with the size of the SSN. The second is a global het- 
erogeneity on the trajectories level indicated by a faster- 
than-linear growth of log2 V with log2 N. Finally, there 
is a heterogeneity on the level of samples of RBNs, indi- 
cated by the fast growth of the cr^ (log2 T)) for RBNs with 
if = 2. 
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